{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Volatility Scenarios\n",
    "\n",
    "Custom random-number generators can be used to implement scenarios where shock follow a particular pattern.  For example, suppose you wanted to find out what would happen if there were 5 days of shocks that were larger than average.  In most circumstances, the shocks in a GARCH model have unit variance.  This could be changed so that the first 5 shocks have variance 4, or twice the standard deviation. \n",
    "\n",
    "Another scenario would be to over sample a specific period for the shocks.  When using the standard bootstrap method (filtered historical simulation) the shocks are drawn using iid sampling from the history.  While this approach is standard and well-grounded, it might be desirable to sample from a specific period.  This can be implemented using a custom random number generator.  This strategy is precisely how the filtered historical simulation is implemented internally, only where the draws are uniformly sampled from the entire history."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, some preliminaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "from arch.univariate import GARCH, ConstantMean, Normal\n",
    "\n",
    "sns.set_style(\"darkgrid\")\n",
    "plt.rc(\"figure\", figsize=(16, 6))\n",
    "plt.rc(\"savefig\", dpi=90)\n",
    "plt.rc(\"font\", family=\"sans-serif\")\n",
    "plt.rc(\"font\", size=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example makes use of returns from the NASDAQ index. The scenario bootstrap will make use of returns in the run-up to and during the Financial Crisis of 2008. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import arch.data.nasdaq\n",
    "\n",
    "data = arch.data.nasdaq.load()\n",
    "nasdaq = data[\"Adj Close\"]\n",
    "print(nasdaq.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the returns are computed and the model is constructed. The model is constructed from the building blocks.  It is a standard model and could have been (almost) equivalently constructed using\n",
    "\n",
    "```python\n",
    "mod = arch_model(rets, mean='constant', p=1, o=1, q=1)\n",
    "```\n",
    "\n",
    "The one advantage of constructing the model using the components is that the NumPy `RandomState` that is used to simulate from the model can be externally set. This allows the generator seed to be easily set and for the state to reset, if needed.\n",
    "\n",
    "**NOTE**: It is always a good idea to scale return by 100 before estimating ARCH-type models. This helps the optimizer converse since the scale of the volatility intercept is much closer to the scale of the other parameters in the model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rets = 100 * nasdaq.pct_change().dropna()\n",
    "\n",
    "# Build components to set the state for the distribution\n",
    "random_state = np.random.RandomState(1)\n",
    "dist = Normal(seed=random_state)\n",
    "volatility = GARCH(1, 1, 1)\n",
    "\n",
    "mod = ConstantMean(rets, volatility=volatility, distribution=dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fitting the model is standard."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = mod.fit(disp=\"off\")\n",
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GJR-GARCH models support analytical forecasts, which is the default.  The forecasts are produced for all of 2017 using the estimated model parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All GARCH specification are complete models in the sense that they specify a distribution. This allows simulations to be produced using the assumptions in the model.  The `forecast` function can be made to produce simulations using the assumed distribution by setting `method='simulation'`.  \n",
    "\n",
    "These forecasts are similar to the analytical forecasts above.  As the number of simulation increases towards $\\infty$, the simulation-based forecasts will converge to the analytical values above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_forecasts = res.forecast(start=\"1-1-2017\", method=\"simulation\", horizon=10)\n",
    "print(sim_forecasts.residual_variance.dropna().head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Custom Random Generators\n",
    "`forecast` supports replacing the generator based on the assumed distribution of residuals in the model with any other generator.  A shock generator should usually produce unit variance shocks.  However, in this example the first 5 shocks generated have variance 2, and the remainder are standard normal. This scenario consists of a period of consistently surprising volatility where the volatility has shifted for some reason.\n",
    "\n",
    "The forecast variances are much larger and grow faster than those from either method previously illustrated. This reflects the increase in volatility in the first 5 days."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "forecasts = res.forecast(start=\"1-1-2017\", horizon=10)\n",
    "print(forecasts.residual_variance.dropna().head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "random_state = np.random.RandomState(1)\n",
    "\n",
    "\n",
    "def scenario_rng(size):\n",
    "    shocks = random_state.standard_normal(size)\n",
    "    shocks[:, :5] *= np.sqrt(2)\n",
    "    return shocks\n",
    "\n",
    "\n",
    "scenario_forecasts = res.forecast(\n",
    "    start=\"1-1-2017\", method=\"simulation\", horizon=10, rng=scenario_rng\n",
    ")\n",
    "print(scenario_forecasts.residual_variance.dropna().head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bootstrap Scenarios\n",
    "\n",
    "`forecast` supports Filtered Historical Simulation (FHS) using `method='bootstrap'`.  This is effectively a simulation method where the simulated shocks are generated using iid sampling from the history of the demeaned and standardized return data.  Custom bootstraps are another application of `rng`.  Here an object is used to hold the shocks.  This object exposes a method (`rng`) the acts like a random number generator, except that it only returns values that were provided in the `shocks` parameter.\n",
    "\n",
    "The internal implementation of the FHS uses a method almost identical to this where `shocks` contain the entire history."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ScenarioBootstrapRNG(object):\n",
    "    def __init__(self, shocks, random_state):\n",
    "        self._shocks = np.asarray(shocks)  # 1d\n",
    "        self._rs = random_state\n",
    "        self.n = shocks.shape[0]\n",
    "\n",
    "    def rng(self, size):\n",
    "        idx = self._rs.randint(0, self.n, size=size)\n",
    "        return self._shocks[idx]\n",
    "\n",
    "\n",
    "random_state = np.random.RandomState(1)\n",
    "std_shocks = res.resid / res.conditional_volatility\n",
    "shocks = std_shocks[\"2008-08-01\":\"2008-11-10\"]\n",
    "scenario_bootstrap = ScenarioBootstrapRNG(shocks, random_state)\n",
    "bs_forecasts = res.forecast(\n",
    "    start=\"1-1-2017\", method=\"simulation\", horizon=10, rng=scenario_bootstrap.rng\n",
    ")\n",
    "print(bs_forecasts.residual_variance.dropna().head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing the differences\n",
    "The final forecast values are used to illustrate how these are different.  The analytical and standard simulation are virtually identical.  The simulated scenario grows rapidly for the first 5 periods and then more slowly. The bootstrap scenario grows quickly and consistently due to the magnitude of the shocks in the financial crisis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.concat(\n",
    "    [\n",
    "        forecasts.residual_variance.iloc[-1],\n",
    "        sim_forecasts.residual_variance.iloc[-1],\n",
    "        scenario_forecasts.residual_variance.iloc[-1],\n",
    "        bs_forecasts.residual_variance.iloc[-1],\n",
    "    ],\n",
    "    axis=1,\n",
    ")\n",
    "df.columns = [\"Analytic\", \"Simulation\", \"Scenario Sim\", \"Bootstrp Scenario\"]\n",
    "# Plot annualized vol\n",
    "subplot = np.sqrt(252 * df).plot(legend=False)\n",
    "legend = subplot.legend(frameon=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subplot = np.sqrt(252 * df).plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing the paths\n",
    "\n",
    "The paths are available on the attribute `simulations`. Plotting the paths shows important differences between the two scenarios beyond the average differences plotted above. Both start at the same point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2)\n",
    "colors = sns.color_palette(\"dark\")\n",
    "# The paths for the final observation\n",
    "sim_paths = sim_forecasts.simulations.residual_variances[-1].T\n",
    "bs_paths = bs_forecasts.simulations.residual_variances[-1].T\n",
    "\n",
    "x = np.arange(1, 11)\n",
    "# Plot the paths and the mean, set the axis to have the same limit\n",
    "axes[0].plot(x, np.sqrt(252 * sim_paths), color=colors[1], alpha=0.05)\n",
    "axes[0].plot(\n",
    "    x, np.sqrt(252 * sim_forecasts.residual_variance.iloc[-1]), color=\"k\", alpha=1\n",
    ")\n",
    "axes[0].set_title(\"Model-based Simulation\")\n",
    "axes[0].set_xticks(np.arange(1, 11))\n",
    "axes[0].set_xlim(1, 10)\n",
    "axes[0].set_ylim(20, 100)\n",
    "\n",
    "axes[1].plot(x, np.sqrt(252 * bs_paths), color=colors[2], alpha=0.05)\n",
    "axes[1].plot(\n",
    "    x, np.sqrt(252 * bs_forecasts.residual_variance.iloc[-1]), color=\"k\", alpha=1\n",
    ")\n",
    "axes[1].set_xticks(np.arange(1, 11))\n",
    "axes[1].set_xlim(1, 10)\n",
    "axes[1].set_ylim(20, 100)\n",
    "title = axes[1].set_title(\"Bootstrap Scenario\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing across the year\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A hedgehog plot is useful for showing the differences between the two forecasting methods across the year, instead of a single day. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analytic = forecasts.residual_variance.dropna()\n",
    "bs = bs_forecasts.residual_variance.dropna()\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "vol = res.conditional_volatility[\"2017-1-1\":\"2019-1-1\"]\n",
    "idx = vol.index\n",
    "ax.plot(np.sqrt(252) * vol, alpha=0.5)\n",
    "colors = sns.color_palette()\n",
    "for i in range(0, len(vol), 22):\n",
    "    a = analytic.iloc[i]\n",
    "    b = bs.iloc[i]\n",
    "    loc = idx.get_loc(a.name)\n",
    "    new_idx = idx[loc + 1 : loc + 11]\n",
    "    a.index = new_idx\n",
    "    b.index = new_idx\n",
    "    ax.plot(np.sqrt(252 * a), color=colors[1])\n",
    "    ax.plot(np.sqrt(252 * b), color=colors[2])\n",
    "labels = [\"Annualized Vol.\", \"Analytic Forecast\", \"Bootstrap Scenario Forecast\"]\n",
    "legend = ax.legend(labels, frameon=False)\n",
    "xlim = ax.set_xlim(vol.index[0], vol.index[-1])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
